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Dedicated to Professor Armin Leutbecher on the occasion of his 80th birthday 

Abstract. This article studies the existence of long-time solutions to the Hamiltonian bound¬ 
ary value problem, and their consistent numerical approximation. Such a boundary value prob¬ 
lem is, for example, common in Molecular Dynamics, where one aims at finding a dynamic 
trajectory that joins a given initial state with a final one, with the evolution being governed 
by classical (Hamiltonian) dynamics. The setting considered here is sufficiently general so that 
long time transition trajectories connecting two configurations can be included, provided the 
total energy E is chosen suitably. In particular, the formulation presented here can be used to 
detect transition paths between two stable basins and thus to prove the existence of long-time 
trajectories. The starting point is the formulation of the equation of motion of classical mechan¬ 
ics in the framework of Jacobi’s principle; a curve shortening procedure inspired by Birkhoff’s 
method is then developed to find geodesic solutions. This approach can be viewed as a string 
method. 


1. Introduction 


The aim of this article is to study the existence and give a consistent approximation procedure 
of the boundary value problem for the conservative dynamical system 


( 1 ) 


dt^ 


-VV(g), 


where V is a smooth potential on Q. We assume that Q is an open subset of K" as this is the 
relevant case for the applications we have in mind; extensions to a more general setting are possible 
but not discussed here. 

For the boundary conditions, we write 


( 2 ) 


9 ( 0 ) = Qa and q{To) = qb 


with qa,qb G Q and Tq > 0. Here, Tq is part of the problem and has to be determined (however, 
the total energy E, defined as the sum of kinetic and potential energy, is fixed). The focus on the 
boundary-value problem is motivated by applications, as discussed below. 


1.1. Hamiltonian systems, rare events and path sampling. Equation Q (furnished with 
various initial or boundary conditions) can be reformulated as the classic Hamiltonian problem 


( 3 ) 


OH 

q = -T^KP^q) 


dq 


for p,q G M", where H is the Hamiltonian 


( 4 ) 


H=-p^ + V{q). 
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Mathematically, the existence of solutions to (§, often more succinctly written as 

(5) ^=(ld 0^) ’ 

with 0 := (p,q), is a classical problem. Periodic solutions have been a particular focus, and 
existence results obtained until the early 1980s are discussed in the beautiful survey article m- 
Already for periodic solutions, a clear distinction has to be made for local results (that is, short 
time solutions) and global solutions describing solutions in the large. Apparently the first global 
global result was obtained by Seifert [22] for a Hamiltonian which is slightly more general than 
the one in Q. The key idea of his proof is based on differential geometry, using an equivalent 
reformulation of in which solutions can be found as a geodesic in a (degenerate) Riemannian 
metric, the so-called Jacobi metric. A curve shortening procedure proposed by G. D. Birkhoff [31 
Section V.7] can then be applied to show the existence of a geodesic. This result has later been 
extended by Weinstein, and a more general result based on a different variational approach was 
given by Rabinowitz [18]. 

In the Sciences, the interest in non-periodic long time solutions has recently been rejuvenated 
by various applications. Namely, complex systems in physics, chemistry or biology can often be 
described by a potential energy landscape with many wells, separated by barriers. A common 
problem is then to find a trajectory joining a given initial point (configuration) with a given final 
point. We study this problem in the situation where the dynamics is determined by Q, and the 
points given in ([^ are potentially far apart. In particular, the two configurations will generically be 
located in different wells of the energy landscape. Rare events are an example of these transitions 
between two wells. Typically, thermally activated reactions have many deep wells separated by 
large energy barriers. Reactants will then spend most of the time jostling around in one well 
before a rare spontaneous fluctuation occurs that lifts the atoms of the reactant over the barrier 
into the next (product) valley. Information on rare events is crucial since they represent important 
changes in the system, such as chemical reactions or conformational modifications of molecules. 
A major challenge in Molecular Dynamics (MD) is that these hopping events take place so rarely 
that the computational limits of MD simulations can be easily exceeded. Since the problem Q-© 
is central in MD, a number of solution strategies have been proposed; see |2()] for a brief review 
of some methods. Further, for practitioners of MD, the question arises whether any numerical 
approximation shadows a physical one |12j (and if so, whether it shadows a generic physical 
trajectory). The lack of hyperbolicity rules out standard tools to prove shadowing (e.g., [TSl 
Theorem 18.1.3]). Thus, for MD, computations are “based on trust” [IHl Section 4.3]. 

One further difficulty for the computation of Hamiltonian trajectories, in particular for MD, 
is that these trajectories are often chaotic, and one has to restrict oneself to averaged statistical 
information. However, of particular interest in the analysis of rare events are trajectories going 
directly from one well of a potential to another one; such transition paths can be used to define so- 
called reaction coordinates. Often the most efficient algorithms using for example path sampling 
assume a knowledge of these coordinates and thus these non-chaotic trajectories are of significant 
practical importance. The method presented in this article is concerned with the calculation of 
such non-chaotic transition trajectories. Even these relatively “simple” trajectories are in practise 
very hard to compute, since they correspond to rare events and take place on very long time scales. 

1.2. Jacobi metric and Birkhoff curve shortening. Rather working with the Hamiltonian 
boundary value problem 0-® directly, we use an equivalent variational formulation, namely 
the Maupertuis principle, according to which trajectories to Q with total energy E are suitably 
re-parametrised geodesics with respect to the Jacobi metric 

(6) g^,^%q):={E-V{q))S,,{q) 

(more generally, if Q is equipped with a Riemannian metric gtj, then gf^^{q) := {E — V{q)) gij{q))- 
So Hamiltonian trajectories are critical points of the length functional associated with (|^. While 
the equivalence has been known for centuries, it seems that little advantage has been taken of the 
fact that this variational formulation has a very convenient mathematical structure. Note that 
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Figure 1. BirkhofF’s algorithm, for the toy example of the Euclidean metric in 
and for i = 5. The initial curve is plotted as a dashed line. Points with odd 
index are marked by black dots, points with even index by grey dots. In the first 
step, the points with even indices are kept fixed, and joined by a geodesic. New 
positions for the points with odd indices on the new curve are determined (solid 
curve). In a next step, these points are joined by geodesics, which determines 
new positions for the points with even indices (dotted curve). The curves (slowly) 
converge to the geodesic line segment connecting qo and q 2 i. 

Hamiltonian problems such as Q are commonly indefinite, while a geodesic problem is elliptic 
and thus bounded from below. 

So the problem reduces to that of hnding geodesics in the Jacobi metric. This problem has been 
addressed by Birkhoff, who described a curve shortening procedure to find global geodesics under 
the assumption that local (sufhciently short) geodesics can be found explicitly. This assumption 
is not met here since local geodesics have to be approximated; the main task is to show that 
nevertheless global convergence can be obtained for a suitably devised local approximation scheme. 
Since the local scheme we propose also relies on a Birkhoff curve shortening idea, we present his 
idea in the global setting first to keep the presentation self-contained. Initially, one joins the given 
initial and final point qa and qt by an arbitrary curve. Then sufficiently many points are marked 
on the curve so that a local geodesic can be computed between next to nearest neighbours (see 
Figure [^. Note that the argument assumes that local geodesics can be computed with sufficient 
accuracy. In the first step, the points with even indices are kept fixed, and joined by a geodesic. 
New positions for the points with odd indices on the new curve are determined (solid curve in 
Figure]^. The procedure of joining next to nearest neighbouring points by local geodesics is then 
repeated for the points with odd indices (dotted curve in Figure[^. It is not hard to show that this 
iterative procedure decreases the length. Under suitable assumptions (e.g., [14]), this method can 
be shown to converge to a global geodesic; however, there are situations such as for degenerating 
metrics where convergence will not take place. 

The central result of this article is an analogous local result, introducing a sequence of approx¬ 
imating sequences converging to a sufficiently short geodesic. The trade-off is that the result is 
local (applicable only to sufficient short geodesics), but does not assume that the approximating 
sequences consist of exact geodesic segments. We rely on the observation that the global Birkhoff 
argument localises the geodesic problem in a geometrically tractable way. That is, one can restrict 
the local analysis to points which are sufficiently close. We then show that they can be joined by 
a geodesic which in addition can be represented as a graph and prove consistency and convergence 
of the proposed local approximation. 

1.3. Results. The main result of this article is a consistent approximation of what we call local 
geodesics; important aspects are that the proof is constructive and yields bounds on the allowed 
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distance between points to be joined by local geodesics. The bounds are not in terms of the 
usual estimates from differential geometry (such as the injectivity radius, that is, the radius for 
which there is a unique geodesic starting at the centre, with arbitrary velocity), but are expressed 
explicitly in terms of the total and potential energy of the molecular Hamiltonian system Q . We 
point out that locality of the geodesics does not necessarily require that the two end points are 
very close; an important aspect of the proof is that we may consider local geodesics which can be 
represented as graphs. Global geodesics can violate this assumption, while suitably small segments 
of geodesics remain geodesics and can be represented as graphs. Birkhoff’s idea to segment an 
original connecting curve allows us to confine our algorithm to such local geodesics. The efficiency 
and applicability of a global Birkhoff method as in |21j will depend on the chosen parametrisation 
(number and location of points in Figure [^. For example, it is possible that refinements in a 
numerical implementation are required. However, the proof shows that in the setting studied in 
this article, no refinement or reparametrisation is required. 

As a by-product, we show the existence of a continuous (physical) trajectory for suitable points 
Qa, qb, using the Jacobi formulation as Seifert [22], but replacing the periodic boundary conditions 
considered there by Equation (i). There are some related existence results [MlE], which also 
rely on the Jacobi formulation. The novelty of the results presented here are twofold: (i) while 
a formulation using the Jacobi metric is natural, a difficulty is that the metric degenerates at 
the boundary dQ of the configuration space, where kinetic energy J dx and potential energy 
/ V (g) dx agree. We provide a priori bounds to ensure that the geodesic stays away from dQ] the 
bounds depend on the location of the boundary points qa and qb or on the total energy E. Bounds 
could be obtained along the lines of thought presented in this paper, but in a simpler fashion. 
The reason why we give a more complicated argument is that the approximation we give here 
is constructive, giving existence of a solution and at the same time a consistent approximation 
procedure. Thus, we obtain an approximation procedure which may not necessarily be the most 
efficient but one for which can dispense with the need for trust. Since the algorithm we develop is 
consistent, the issue of shadowing is answered in an affirmative way for the procedure we propose, 
under the assumptions made on the potential energy V and the bounds in term of total energy E 
made on the end points. 

The existence of geodesics joining a given initial and final point in open domains Q, which 
is trivial within the radius of injectivity, is not obvious if the two points are further away from 
each other. We point out that the argument of this paper automatically proves the existence 
of such a geodesics for the Jacobi metric with sufficiently large total energy E. Other existence 
results, proceeding along quite different lines, can be found elsewhere mm- We remark that 
we will not address the question of how to choose E] this choice typically requires insight in the 
physics, chemistry or biology of the problem in question and thus cannot be answered in the 
general mathematical framework considered here. However, the existence result given here can be 
interpreted in two ways: given the initial point qa and the total energy E, the arguments provide 
estimates on possible locations of the final point qb so that qa and can be joined by a trajectory 
with total energy E. Alternatively, given g^ and qb, the analysis provides lower bounds on E such 
that qa and qb can be joined with this total energy. It is easy to see that no general existence 
theorem can hold if qa, qb and E are unrestricted {E determines the configuration manifold Q, 
and in particular for small E the configuration manifold may be disconnected, so qa and qb could 
be in different components). 

1.4. Applications and limitations. Approximations which are proven to be consistent, such 
as Godunov’s scheme for hyperbolic equations, are often less efficient than algorithms for which 
consistency cannot be shown. The approximation introduced in this article is no exception. Yet, it 
is often possible to take inspiration from a consistent approximation and deduce efficient (but not 
provably consistent) formulations. This is the case for the formulation introduced here; a related 
flow model approximation has been shown to be able to detect different trajectories joining points 
in different wells of the energy landscape of the Muller potential and the collinear reaction H 2 -I-H — 
H -b H 2 |20|. While our formulation relies both on the choice of Maupertuis’ formulation and 
Birkhoff’s curve shortening, which seems to be new to the field of numerical methods for Molecular 
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Dynamics, the curve-shortening procedures resembles other rubber-band algorithms m- For 
the isothermal case, a nudged elastic band method has been proposed by E, Ren and Vanden- 
Eijnden [^. There, the aim is to find minimal energy paths, which are defined as paths along 
which the orthogonal component of the deterministic vector field vanishes. The approach in to 
reduce the orthogonal contributions iteratively bears many similarities with the method presented 
here. As examples of string methods with temperature, we refer to the pioneering work by E, Ren 
and Vanden-Eijnden mm- 

Maupertuis’ principle has been used before in MD [T], but without the connection to Birkhoff’s 
curve shortening method. We also refer the reader to the recent work by Cameron, Kohn and 
Vanden-Eijnden which gives an analysis and in particular convergence results for a steepest descent 
string method 

The usual numerical approach for a boundary value problem is a shooting method; there, the 
existence of a solution is assumed as well as the closeness of an initial guess to the solution. 
(There are abstract existence results available, see for example [TOj; however, as noted by Stoer 
and Bulirsch [231 7.3.3], the abstract formulation of the boundary conditions to be imposed rules 
out the condition ii). for the first-order system considered there, even for the case n = 2.) The 
method discussed in this paper provides both an explicit existence proof and estimates on the 
closeness required. 

We point out that the Birkhoff algorithm also provides a strategy for gluing together local 
geodesics to obtain global paths, which converge to a global geodesic. Note that this latter aspect 
of the Birkhoff approach defines a natural tool to localise the computation of geodesics which can 
be exploited for parallelisation. This aspect is discussed in more detail in |5T]. An advantage of 
Birkhoff’s method is that these local steps are intrinsically parallelisable. 

Einally, we mention connections to the Onsager-Machlup / Freidlin-Wentzell theory. There, an 
action functional is derived, with the minimal action path describing the most likely trajectory. 
The connection between that theory, Hamilton-Jacobi theory and the Maupertuis principle as 
discussed here is a topic with many open questions; we refer to |5] for results in this direction 
together with applications to rare events. 

1.5. Notation. Throughout the presentation, Q is the configuration manifold of a system and 
thus describes all possible states the system can occupy. The coordinates of the phase space 
(cotangent bundle) T*Q are (qfipj), position and momentum. Analogously, the coordinates of 
the tangent bundle TQ are where q^ denotes the velocity. We assume that the system 

dynamics is conservative with 3V degrees of freedom. Then, the Hamiltonian H: T*Q ^ M. is 
defined as 77 := E := T + V. Here, the kinetic energy T = T{p) is a function of the momenta only 
and V = V{q) is the potential energy, depending on the coordinates q alone. The Lagrangian of 
the system is a function L: TQ —)■ K, namely L{q,q) = T — V. For a wide class of applications, 
it is sufficient to consider {q,q) = inner product for a system with N particles 

with mass mj , and to assume that the Lagrangian is of the standard form 

(7) L{q,q) ■=^{q,q)-V{q). 

This paper is organised as follows: In Section we review the Maupertuis principle and the 
Birkhoff curve shortening algorithm, both for the continuous setting. Section describes the 
analogous discrete setting, contains the relevant a priori estimates and introduces the discrete 
Birkhoff procedure for a fixed discretisation. Section [^ describes the Birkhoff refinement and 
convergence to the continuous limit. Numerical simulations and numerical convergence rates for 
a model problem are the content of Section]^ 

2. The continuous setting 

Our construction will be guided by a variational formulation, equivalent to Q, where convergent 
approximations can be obtained with relative ease. This continuous setting is sketched in the 
present section. 
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It is a well known fact that solutions to 0 with pre-assigned total energy E are re-parametrised 
geodesics in Jacobi’s metric 0 . This formulation is sometimes denoted Maupertuis ’ principle, or 
Jacobi’s least action principle. For the special Lagrangian Q, the Routhian associated with 
Jacobi’s principle is 

( 8 ) R{q,q') = 2{E-V{q)){q',q'). 


Obviously, is a metric in those regions of Q where V(q) < E. The action functional 

pb 

(9) 


J[y] ■= f i?(7,7')dr 

J a 


is the measure of the length of 7 in this metric. For a given curve 7 , the value J[ 7 ] is often called 
the energy of 7 ; the length of the curve is then 

( 10 ) 


L[l] ■■= f VR ( 7 , 7 ') dr. 

J a 


It is trivial to verify that critical points of the energy functional are critical points of the length 
functional; the converse is true for curves parameterised by arc length. The length functional is 
invariant under re-parametrisations, while a minimiser of the energy functional is automatically 
parameterised by arc length. 

The Maupertuis’ principle seeks geodesics, that is, stationary solutions of the functional 0 
with the metric ([^. Maupertuis’ principle has been employed in a number of computational 
approaches and is regarded as a very accurate method for the verification of other algorithmic 
formulations m- 

We already mentioned in passing that solutions of 0 are re-parametrisations of solutions of ® 
(or, equivalently, @). The re-parametrisation is such that the physical time for a solution of Q 
can be recovered via the explicit formula 


( 11 ) 


t = 



2{E-V) 


ds. 


With the exception of [20], numerical methods for Q seem, to the best of our knowledge, not 
have taken advantage of the geodesic formulation (|^ (see [ 20 | for a recent survey and a method 
that relies on observations similar to, but simpler than, those made in the present article). This 
is somewhat surprising, since the Birkhoff curve shortening algorithm is a classic method for the 
convergent approximation of geodesics. 


2.1. Existence of extended geodesics. Birkhoff’s curve shortening method |3] is a construc¬ 
tive way to find extended geodesics, based on the assumption that local (short) geodesics (that is, 
geodesics joining points within the radius of injectivity) can be computed exactly. In this subsec¬ 
tion, we recall the classic Birkhoff method; the main part of the paper then addresses the question 
of how to find the local geodesics constructively, for the case of the Jacobi metric. A straightfor¬ 
ward implementation of the Birkhoff method relies on an approximation of local geodesics within 
the radius of injectivity and thus requires knowledge of the radius of injectivity. For the complex 
energy landscapes we have in mind, this radius is not easily computable in a quantitative way. 
Thus, we have here two aims. Firstly, we present an algorithm for the numerical approximation 
of local geodesics in an explicitly given neighbourhood. Secondly, we obtain a quantitative de¬ 
scription of the size of this neighbourhood. A Birkhoff method then glues together these local 
geodesics to obtain an extended piecewise geodesic curve. 

The aim of this article is to develop a discrete framework that mimics the Birkhoff procedure. 
Besides the usual difficulty of discretisation errors inherent in any numerical approach, we face 
the challenge that even the computation of local geodesics, as in Birkhoff’s algorithm, is time- 
consuming and difficult to control for non-Euclidean metrics. We propose an approximation of 
Jacobi’s metric Q by the trapezoidal rule. The key observation is that the difference between 
the Jacobi metric and the Euclidean metric occurs on a fine scale, described in greater detail 
below. The analysis of Section will show that crucial bounds on the discrete curvature in 
a Birkhoff procedure can be obtained since it is possible to show that in some (quantitatively 
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characterised) situations the Birkhoff procedure for the Jacobi metric is locally identical to that 
for the Euclidean metric. The analysis will also reveal that for other configurations, that is, other 
geometric configurations, the two approximations differ, which results in the Jacobi procedure 
making steps which seem counter-intuitive if regarded within a Euclidean picture. Obviously, 
such a disagreement of the two approximations is necessary as we need to compute a geodesic in 
the Jacobi metric and thus have to differ at some point from the Euclidean picture. 


3. A LOCAL DISCRETISED BiRKHOFF METHOD 
This section mimics the continuous framework laid out in Section in a discrete setting. 


3.1. The discrete setting. Throughout this section, we assume that the total energy E is suf¬ 
ficiently large, as described in the next paragraph. We point out that the choice of E determines 
the configuration manifold Q C K", which we take as the set of points q where E — V{q) > 0. 
Let qa and qb £ Q he given; define £ = , where I-] is the Euclidean distance on Q (not the 

Jacobi distance). We choose an orthonormal basis for Q such that qt — qa = 2£ei. For q £ Q, we 
write q = (A, T) G K x and in particular ei = (1, 0). 

Let us write the Jacobi metric in the form 


( 12 ) 


gfriq) = 


We require E to be sufficiently large so that the line segment joining qa and qb is contained in 
Q as well. In fact, we will work in a framework where either E is chosen large enough, depending 
on the given points qa and qb, or, given E, we choose points qa and qb in Q with sufficiently small 
distance £, such that qb £ B({qa) C Q. 

We consider a convex set Qh such that Qh is compactly contained in Q, Qh (H Q', then the 
Jacobi metric is not degenerate, and hence Riemannian, on Qh- It is bounded from above and 
from below by Euclidean metrics, but we do not use this fact. Then, if V is sufficiently smooth, 
the metric factor h is C^’'^{Qh) and there is a finite iL G K such that the estimates 


(13) 

(14) 


<iL(|x| + H)'+ffi 

VnKX + x,Y + y)- VnHX, Y)\ < H (| x | + \y\f 


h{X + x,Y + y)- { h{X,Y) + x — h{X,Y)+yXNHX,Y) 


hold. 


Applying Equation (10) to the length of a straight line segment [ 91 ,^ 2 ] between two points, 

7[9i.92](0 ■=qi+t ((?2 - qi) for t £ (0,1), 


we obtain 

(15) = / e'‘(7[9i.g2](^))||7[gi.,2l(^)|| di = / e'‘(7[g^,,,](t))dt Ijgi - g2||- 

Jo Jo 

We now introduce the discretised setting. We first define an equidistant Cartesian grid on 
Qh C M". We discretise the integral for the length of a straight line segment as follows. 


Definition 3.1 (Discretised length of segment). For the straight line segment {qi,q 2 ] between two 
points on the grid, we define the discretised length by applying the 2-point trapezoidal rule to (15) 


■^[ 91 . 92 ] •" 


e^(7[9i,92](0)) +e^(7[9i,92](l))| 


qi - <72|| = 


e^{qi) + e^{q2). 


\qi - 92] 


Below we will introduce a discrete Birkhoff procedure which chooses to move points of polygonal 
curves in a direction normal to ei in order to decrease the length. It is thus necessary to estimate 
changes in the length as normal variations of a curve are considered. This is the content of the 
following lemma. 


Lemma 3.2. Let Qh be convex set such that Qh <21 Q and assume the metric factor satisfies the 
Holder estimates (13) and (14) for H > 0. Let q = {X,Y) £ Qh be a grid point. 
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Assume further that e 0 and S, A G M" ^ are sueh that <5 7 ^ 0 and ?£ = 9 + (s, A) and 
qs = q + (0, S) are also grid points in Qh satisfying the bounds 


(16) 

Then the difference quotient 


A 


6 


< 1 and 


£ 


£ 


(17) 

satisfies 


D{e,S,A) := 


< 1 . 


^[q,gs] 


|e<5| 


^ f e\q.)F{^) ^ 1, 

l<J| 


D{e,5,A) = ^-\-AAlAdA^^lS7Me^{q)PF[-] 1+0 


■Oder), 


where 


(18) gf = = g + l(e, A) and Ppir]) := + ^( 77 ) := VPf(??) = 


Pf{v)' 


Note that |+’( 77 )| < 1 for all 77 G 
Proof. We can write 




[95,<7e] 


e^jqs)+ e'^{qe) ( A - 5 


Then it follows directly that 


D{e,5,A)=DiS2 + SiD2 


with 


r e^(g^) + e'»(gr e'^(g) + e^(g,) 
2 2 
e'*(gj) + e'‘(g£) , e^{q) + e^{qe) 


Si = ; 


= ^ K(g5)-e^g)] , 

= 2 [^^( 9 ) + e^(ge)] + ^ ■ Oi, 


= M (^) - (7 

= I ( 7 )) = ( 7 ) + f 


Let us define Diq by writing 


Oi = -|^ • VNe^iq) + Dio, 


where we estimate the error term Diq, using the Holder bounds for the metric factor h 

\Dio\ = ^|e"(gj) - e\q) - S ■ V^ve^g)] < \Sr^. 


Secondly, we consider Si and write for the first term 


1 


[e\q) + e\q,)] = e\q.J + S 


10 


with qc defined in (181 being the mid-point on the segment [g, g^]. By symmetry, the Holder 
bounds imply 

|5io|<H(|e| + |A|)i+“=i7 (1+ '' ^ 


|ed+“ < 2H|ed+“. 


Furthermore, we deduce 
1 
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where F is defined in (18). Thus, we can write 


= — ITT • -P + D20, 
|(5| e 


where, using \DF(ri)\ < C for all 77 S M, 


-D20I = 


1^1 


0 JQ 


£ 


DF 


A — 6t 


£ 


dr ds 


< 


C. 


We summarise 


D1S2 + S1D2 = ( ^^ • Vve^g) + + ^^2 

+ ( e^{qi) + 5 io + f Pi) • \f (^^) + P20 




1 . 


fx\\ ^ 

5 

- 


J J 



+ o(N“), 


with the functions F and Pp from (18). 


□ 


We now prepare a crucial quantitative upper bound for the discrete bend of a polygon (Propo¬ 
sition 3.3). To this aim, consider three neighbouring points along the polygon. We assume that 


for this triplet (g_e,g,g+£), the X co-ordinates are distributed equidistantly for an £ > 0, 

q={X,Y), = q +{±e,A±). 

For a given 5 7 ^ 0, we want to estimate the centred differences 

AL{e,d) := (T[ 5 _j,g+(o, 5 )] + P[(j+(o, 5 ),( 3 +e]) ~ {L[q_^^q] + £[ 5 , 9 +^]) • 

Using centred coordinates, we can rewrite this as 

AL{e,S)=e\S\ (5" + P+) , 


where 


:= P(±£, (5, A±) 


with Zl(±£, 5, A±) defined in (17). 

We will now combine the length calculation of Lemma |3.2| on both sides of the centre point q. 
We define, in analogy to (IT^, 


(19) 

( 20 ) 


g±£ = g + |(±u A±) 



±e: 
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and obtain 

AL{e,d) 


= D- +D+ = 

S I I 

M'l-^ + + 


—e 


( 21 ) 


+ 0 
6_ 


+^V^e^g)Jl + 


o(kr) 


A, 


e^{q^)F+-e^iq^)F_ ^ A+ ^ 

-a^ + Vjve'‘(g)Wl+ + 




1 + 


A4 


— 1 / 1 + 


A 


V 


+ 0 


+ 0{\en 


M 


eHqi)F+-eHq^)F- 


+ ^ NG^{q)\ 1 + 


A+ , A_ \ / A+ A 


l<5| 


O 


Let us remark that the term 


^Ne\q)- 


6 —e / \ € —£ 


2 1 1 


o(kn. 


1 + 


A 


is a proper difference quotient for the discretisation length e. 

The next result shows that the length of the edges of a three point polygon can be reduced by 
moving the middle point towards the line connecting the two outer points. 


Proposition 3.3. Let ^ be the l°°-sphere in 


pn—1 


:= <! ^ e 


pn—1 


sup \ej • < 1 and \ek • = 1 for some fc S {2 ,..., n} > . 

jG{2,...,n} I 


There exists N > 1 and eo = Sq{N) S (0,1/A) such that for all e G (0,£o) o.nd all triplets 

{q-e,q,qe) with q±e = q +{±£,A±) 

which satisfy 


( 22 ) 

and 


Ah 


±£ 


< 1 


A+ + A 


(23) 


Nv, with N G (A, 3A) and v G see Piff. 
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Figure 2. Geometric configuration as set out in Proposition |3.3| 


there holds 

j AL(e, (5) > 0 for every S with |5| < e“+^ and ||| • = 1, 

|AZ(e, (5) < 0 for every <5 with |5| < e“+^ and ^ ■ v = —1- 

The important implication of this statement is that length-reducing procedures for triplets will 
not increase the discrete curvature indefinitely. Specifically, if the curvature of a triplet is such 
that it is in the white inner square in Fig. then the length shortening procedure may increase 
the curvature (unlike in the Euclidean case). However, if the curvature increases, it will eventually 
enter the grey region in Fig. Then the proposition shows that a further step “outwards” (in¬ 
creasing the discrete curvature) necessarily increases the length, while the corresponding “inward” 
step decreases the length. This prevents the discrete curvature to grow without bounds under 
a length shortening process. So in the inner white region, the Birkhoff procedure for the Jacobi 
metric and for the Euclidean metric can differ; in fact they have to differ at some point since the 
results, the respective geodesics, differ. 


Proof. From (21) we deduce, as |^| < |e|“, 

e^{qe)F+-e'^{q^)F_ 


AL{e,S) S 


+ ^ N^^{q)\ 1 + 


(24) 


We rewrite 


l<5| 


V^e^g)- 


A+ , A_ \ / A+ A 


6 —£ / \ € —e 


2 if 


0{\er). 


+ \ 1 + 


A 


, e^(9^) + e^(g^) F^F 

(25) e\qs)F+ - e\q^)F. = -^ (F+ - F_) + " 


(e\q^)-e\q^)). 


Firstly, with the identities 
(26) 


A± _ A+ - A_ ^ A+-b A, 
±e 2e 2e 
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we infer for the difference of F± 

A+ A+-A_ I A++A_ A+-A_ _ A++A_ 


(F+-F_) = 


A+ +A_ 
2 e 


2e 2e 


2e 


2e 


1 + 


1 + 


A 


1 + 


A4 


1 + 


A 


1 + 


1 + 


A 

—e 


A+ - A_ 

- 1 - ^ 



1 


\ 

2e 

/ 

v\/^ + 

A+ 

£ 


A 

— £ 

2 

/ 


A++A 


A+ 

e 

+\ji+ 

A_ 

— e 

2 

A+-A 

1 |2 1 A_ |2 

1 £ 1 1 — £ 1 

£ 



2 


“ ( 



where 


(27) 


A = 



= A 


A+ + A 




1 + 


1 + 


A 


/ 

A 4. 

2 / 

A 

2 


£ 

+G+ 

— £ 



Secondly, (13), (14) and ( [^ imply 

e'^(q) 


<ei+“iJ-2i+“+e||V/i|| -2. 


These two steps imply for (251 


e'^(qs)F+ - e\q^)F. --^ 71 ^++ ~ 


< e'^iq) {2^+°‘e^+°‘H + 2e\\Wh\\) 

V 2 


where we used the bound |F^| < ^ implied by ( 22 ). 


We return to (24), and with A from (27), together with (23) and 
we obtain 


A+-A 


At + 

£ ' —£ 


< 2 , 


AL{e,6) ^ e''(g|) + e'*(g^) 
2 


(5 A++A_ 


-livhii 


2^2 


- (l + |jVh||.2^) -||Vh||y2 

+ o(N“) 


A+ + A 


> 


eHqi) + ^Hq^) 


2 \ 


—Ap • iV - 1 + ||Vh|j • ^ - \\Vh\\V2 


1^1 

-eA||V/i|| 


V2j 


V2 


We use the inequality 


\/l + |a|^ \/l + i^r 


+o(kr 


> 1 Uibi 


for all a,b G ^ to deduce for a diagonal element of the matrix A 


efAe, > 


1 + \ I + 


1 + 


A„ 
— £ 


- A 


72 " 
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Thus, we obtain in the case ||| 


AL{e,6) 


i) = l 

e'‘(94)+e'*(9^) 


> 


Hence 


e|5| 

AL{s,S) 

e|5| 


V2^ 

e'*(9e)+e'‘(g^) 


{iV - \yf + ||Vh|| • 2 (2 + 2 + eiv)] } + Oder). 


> 


(tv - 2 (4 + ||V/i|| [4 + iVe])) + 0(|e|“). 


Thus, we may take N = 3 (4 + 5|| V/i||) > 12 and choose eo smaller than 1/iV. As N > N and as 
we may reduce eo possibly further to compensate the error term, we can ensure the positivity of 
for all 0 < e < eo- 


Analogously, we deduce in the case ||| ■{> = — I that there holds 


AL{e,6) 




< 


eri - 2 

that is, strict negativity for the opposite sign. 


(-7V + 2(4+|lVh|| 


4 + TVe 


+ O(|er)<0, 


□ 


3.2. BirkhofF method for a fixed discretisation. We recall that, for given qa and qt and 
£ _ l9°~9i’l ^ where dl is the Euclidean distance on Q (not the Jacobi distance), we have chosen an 
orthonormal basis for Q C K" such that qb — qa = 2f'ei. For q £ Q, we write q = (A, Y) £ 

Definition 3.4 (Polygon, associated points and differences). Let M £ N large enough such that 
s := ^ < e. We define a polygon with 2M + 1 vertices as 


= -M,..., M. 


1 / \ 

(28) qj = {Xj,Yj) :=-{qb - qa) + iXjei+ '^Ykek+i\ with Xj = je for j 

Note that qa = q-M = (A_m,0) = {—£,0) and qb = qn = (Am,0) = (£, 0). The polygon 7 
associated with these points consists of the line segments joining neighbouring points. 

We then define the set of interior nodes ff := {j £ Z \ \j\ < M} and set for all j £ J 


Remark 3.5. Note that 

A" 


Af = - y,. 

and 

A+ = y,_i + y,+i-2y,. 


We now show that the fixed boundary points q±M and the estimate on the second differences 
from Proposition |3.3| ensure that the first differences remain bounded by 1. 


Definition 3.6. Consider a polygon 7 as in Definition \3.4\ Let N be given by Proposition \3. 
We say that 7 G */ the second difference quotients of 7 satisfy 


(29) 

Lemma 3.7. Let t < Iq '.= 


max max 

j^J i=2,...,n 


A: 


a; 


< 2N. 


4Afv 


where N is given by Proposition 3.3 Let j be a polygon as 


in Definition 3.4 If J £ D 2 N (see Definition\3.0^, then there holds 


(30) 


sup 


A^ 


±e 


< 1 


as well as 7 is contained in 

(31) r2N ■■= {g = (A, r) : A G H, i] ,\Y\<N - jAj^) } . 
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Proof. We add up the second differences: As 1 ±m = 0, any unit vector e^, /c = 2,..., n, orthogonal 
to ei satisfies Ck ■ Y±m =0. Hence, for A+ = Yj+i — Yj it follows that 

0 = ek ■ {Ym — Y-m) = ^ efc-A+. 

By the Mean Value Theorem, there exists a j- G J U {—M} such that ej, • A^_ < 0. For any 
j G J U {—M} we thus deduce from the boundedness of the second difference quotients 

Cfc • A+ = Cfe • A+_ + sign (j - j_) ^ Ck ■ (A^^ - A+) 

/=min{j,j_} 

= Cfe • A+_+sign(j - j_) efe-(A+i + A[;J 

Z=min{yj_} 

< \j -j-\- 2Ne^ < 2M ■ 2Ne^ = 2t ■ 2Ne. 

We proceed analogously for a j+ G U {—M} with • A^ > 0 and obtain for all j G JU {—M} 




< 


As Cfc is arbitrary, we find 


< AN£y/n — 1. As £ < £o, we conclude 




< 1 . 


Note 

□ 


As A^ = — Aj_|_;^, the estimate extends to the corresponding ‘negative’ differences A^ 
that an integration of the second differences yields directly the inclusion 7 G V^n- 

We now define the Birkhoff method for fixed e > 0. We start with a polygon 7 ° represented 
by the vertices that can be written as in ( |2^ . We think of the Birkhoff method as an iterative 
process to update a polygon 7 * to a polygon 7 *+^ which has a strictly smaller discrete length. It 
will be shown later that such an update is, for a fixed discretisation, not always possible and the 
Birkhoff method thus terminates. 

Definition 3.8 (Birkhoff step). We consider a polygon 7 ^ represented by 2M + 1 points gj, as in 


Definition 3.4 Let 0 < ^ < and 

N := {aci I a G {±1} and i = 2,... ,n}. 

Then 

(1) we consider sequentially every j G J. For given j, consider sequentially every v G N and 
try to move the interior point gj to q* = qj + 6, with S := C,v, to achieve 

That is, the passage via q* is shorter than via the original gj. In the affirmative case, then 
we define the update 7 *+^ as 

9 ^+^ := q* and := g[ for k yf j. 

Thus the update has strictly smaller discrete length. 

(2) If \(1)\ is not affirmative for any j G J, then Birkhoff step is called void. 

The Birkhoff step depends on the sequential order chosen for J and Af. Here and later, we regard 
this choice as fixed and hence the Birkhoff step is uniquely defined. 


It is immediate that if g) is on the grid defined in Section 3.1 then lies on the grid as 


well; the Birkhoff step thus makes only movements which are compatible with the grid, and thus 
results in polygons with vertices on the grid. 
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Definition 3.9 (BirkhofF method and map). The Birkhoff method is the iteration obtained by 
consecutive Birkhoff steps starting with a polygon 7 ° until the Birkhoff step is void. 

The Birkhoff map maps a starting polygon 7 ° to the final polygon obtained by the Birkhoff 
method. 


Proposition 3.10. Let £0 > 0 be given by Lemma \3. 7| with N > 1 be given by Proposition \3.3[ 
Consider an initial polygon 7 ° as in Definition 3.4 with sufficiently close endpoints, that is £ < Iq. 
Furthermore, assume 7 *^ G D 2 N, that is, the second differences are bounded. The updates of the 
Birkhoff steps obey the same assumptions, that is the remain graph-like polygons in the sense 
of Definition \3.4\ and are contained in D 2 N ■ Furthermore, the Birkhoff method terminates after 
finitely many steps with a final polygon in the strictly smaller set Dn > that is, there holds (recall 
Definition's^ 


max max 

i=2,...,r, 


A; 


e* 


a; 


< N. 


We point out and will use later that the second differences bound associated with Djsi is exactly 
half the bound associated with D 2 N- This will be crucial to compensate for the doubling of the 
discrete curvature if one halves the stepsize. 


Proof. Fir stly, we notice that 7 ° G 'P 2 N by Lemma 3.7 Consider a Birkhoff step according to 
Definition |3.8l (i) If the step is void, then 7 ° is the final polygon. The final polygon has to lie in 


I?AT as otherwise Proposition |3.3| would ensure the existence of a further affirmative Birkhoff step, 
that is, an update site qt and an associated direction 6 such that AL{e,6) < 0 . 

(ii) Otherwise, consider a single affirmative Birkhoff step acting on a polygon 7 * G D 2 N, that 
is, there is & j G J and v G N such that 




< u 


For q := qt, q^ := and 5 = (v, it follows that 

(32) AL{e,S) = (T|q_ + T[q+( 0 ,( 5 ),Q'+]) ~ q,_|_]) < 0 . 

We want to show that the Birkhoff step leaves D 2 N invariant. 


Case 1: If 


max max 






< 2iV, 


then by definition of D 2 N the update is contained in D 2 N, so nothing is to be shown. 

Case 2: Now we assume on the contrary that is not in the set D 2 N- That is, there exist a 
j G J and an f G {2,..., n} such that 


A 


— ,/ + l 


A+,z-ti 


thus we can write 




,/ + l 


> 2N-, 


S - C 

-A-2^ =: -Nv-2^v. 


In view of the notation used in Proposition |3.3[ we write 


where A± := A^’* with 


A± 


±e 


< I by Lemma 


3.7 


Furthermore, there holds N G {N, 2N], v G S[ 


n—1 

00 


and S = (v, v G N. Note that v ■ v is restricted to the values —1,0,1 by the definitions of S" 


and N. By assumption, however, 7 ^+^ has left the set D 2 N, hence v ■ v = 1. Proposition 3.3 then 


implies AL{£,S) > 0. This contradicts (32), so Case 2 is in fact impossible. 

Using the inclusion D 2 N C S, we conclude that for a finite discretisation length s, the number 
of distinguished polygons represented on the discrete grid is finite. Each affirmative Birkhoff step 
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is strictly reducing the length, thus does not allow us to visit the same polygon twice. Hence the 
method terminates after finitely many affirmative Birkhoff steps. As already argued in (i), the 
final polygon has to lie in as otherwise one could show the existence of a further affirmative 
Birkhoff step. □ 


4. Birkhoff refinement 


This section consists of three parts. In the first part, we define a sequence of polygons 7 ^,, which 
are the the final polygons of the Birkhoff method for the discretisation length e^. As Sk —>■ 0, we 
will show that the curves 7 ^, converge to a curve 7 . In the second part, we recall a weak formulation 
of the geodesic equation. In the third part, we show that the limit 7 satisfies this weak geodesic 
equation, hence 7 is smooth and a stationary curve for the Jacobi length functional. 

4.1. Refinement and convergence. 

Definition 4.1 (Refinement). Let ALq be large enough such that Eq := ^ < i and Co := Eq. Let 
us define the starting polygon 70 as the Eq- discretisation of the straight line segment with endpoints 
Qa, qt- 

For k = 1,2,..we want to halve the discretisation length, that is, e^ '■= for := . 

Let Jk '■= {j € Z I IjI < Mk\. We embed the polygon jk-i into the finer grid by introducing new 
vertices at the midpoints of the connecting line segments. In the notation of Definition \3.4[ this 
means the embedded polygon fik has the vertices 

for even j e JkU {±Mk}, 

= 5 for odd j G Jk- 

We define Cfc := |Cfc-i- ^nr this choice, the vertices of jk He on the finer grid (e^Z, • 

Then 7 ^ is the Birkhoff map of jk for the discretisation length Ek ■ 

To simplify the notation, we write on each level k 

AJ := Aj~, which implies — AJ_;^ = A~. 

We now want to prove that the Birkhoff refinement laid out above will lead to a converging 
sequence of polygons if we start the iteration with properly chosen initial points. Specifically, 
let us first assume that the total energy E, and thus Q, is fixed, and let us choose a nonempty 

, N = 3(4 + 5||Ve^||) with 

the norm taken over the set P. 


compact set P C Q. Then let N be as in the proof of Proposition 


3.3 


Definition 4.2. Any pair of endpoints {qa, qb) is admissible if they meet the following two condi¬ 
tions, which depend on N: 


(1) i = < £q, with i(j := ]—- as in Lemma 

( 2 ) the set D 2 N defined in Lemma\3. ?| is contained in 



Theorem 4.3. For an admissible pair {qa,qb), we consider the sequence of graphs {x,fk{x)) with 
X G {—i,€) representing the polygons 7 ^ obtained by the Birkhoff refinement. As 7 ^ G Dn, the 
centred differences satisfy for all k > 0 the combined estimates 


(33) 

(34) 


max 

j£JkU{-Mk} 



A'f - Al 

max max • ■—-— 

jEJki=‘2,...,n Ef, 


< N. 


Proof. Observe that the initial polygon 70 is contained in since all its finite differences vanish. 

By induction let k = 1,2,.... According to Definition 4.1 we embed 7 fc-i into the finer grid 
of size Ek, which is half the size of Ek-i. The recursive definition of Cfc implies Cfe < ^'k^°‘^ with 
«=D||-2e(0,i). 
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With respect to the embedded polygons ')k-i are in I?2Ar; this follows since all newly 
introduced vertices have a vanishing second difference quotient —— 2 at odd nodes j G 

whereas all quotients with even j G J'k are doubled in size, 


1 

J j-i 


A 


= 2 - 


fc -1 

‘i-1 


' fc -1 


From Proposition 3 . 10 [ it follows that the final polygon 7^ obtained from the Birkhoff map is in 
fact in the smaller set 'Dn- This is crucial for our argument, since it shows that one application 
of the Birkhoff map after a refinement retains the same discrete curvature bound as before. 

This shows that every polygon 7^, with step-size belongs to Vf^. By Lemma [XT] it obeys also 
a bound on the first difference quotients, which is uniform in k. □ 


Corollary 4.4. For an admissible pair {qa,(lb), we consider a sequence of graphs {x, fk{x)) with 
X G (—£, €) representing the polygons 7^, obtained by the Birkhoff refinement. Then the functions 
fk converge in ((—£,£);to a limit f G ((—£,£);. We write 7 = {x,f{x)) for 
the limit graph. 

Proof. As "fk G 1 ]; 5 ), we deduce the claimed convergence in ^^([ 0 , 1 ]; S') for any fi G ( 0 , 1 ) 

by Arzela-Ascoli. □ 


Let us remark that the graphs of all polygons 7^ and hence the limit 7 belong to the subset 
'P2N of Q defined in Lemma [ 3 . 7 | 


4 . 2 . Variational formulation for a geodesic. Let us consider a geodesic which can be repre¬ 
sented as a graph. In this subsection, we derive a weak formulation for its governing equation. 
Thus, we consider a curve represented as a graph 7(0;) = (cc, f{x)), x G {—£, £). The length of this 
curve in the Jacobi metric yff^iq) = e’^^{q) 6 ij is given by 


For a orthogonal perturbation 7^ = (x, f + ey), where y has compact support in {—£,£), we deduce 
via integration by parts 


d 

de 


Lbe] = f 

2/(x) • + 1 d dJb) ■ dxVb) 

£=0 d-i 

-1 

'I 

-f 

1—1 


dx. 


Thus, a graph 7(x) = (x,/(x)) with f G {—£A) is stationary for the length functional L..^ if 
the variational derivative vanishes for all functions y G Hq (—£,£). We will use this 

weak formulation, rather than the more common one obtained by a further integration by parts. 


4 . 3 . Characterisation of the limit of the Birkhoff refinement. Using a smoother interpo¬ 
lation of the grid points representing the polygons of the Birkhoff refinement we obtain a better 
convergence and a smoother characterisation of its limit than the previous result of Corollary | 4 . 4 [ 

Theorem 4.5. For an admissible pair {qa,qb)> consider the sequence of polygons "fk = (x,/^(x)) 
obtained by the Birkhoff refinement. The functions fk converge in F[^ (—£,£) to a limit f. Fur¬ 
thermore, f G {—£,£) o-nd the limit graph "f(x) = (x,/(x)) satisfies 


( 35 ) 


y{x) ■ -k 1 1^ f{x)f + dxfbb^x^ 



A + Ie/wIJ 


for every function y G JLq {—£,£). 
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Proof. Analogously to the definition of AL(£, S) in (21 1 , we define 

— “f Qj + l\) — 'P T Qj + l\') * 

It is convenient to write for j G Jk in analogy to ( |T^ 

.v‘ 4 :=Av‘ + .Y‘„). 


In further analogy to (20 1 , let 
(36) 


F^:= 


_ j_ 

Sk 


! 


2 

\ J + 



V 

£k 



We want to examine AL{e,d) of (21) evaluated at a node qj for e = and <5 = where 

a G {±1} and i G {2,...,n}. Theorem 4.3 provides the necessary estimates for the difference 
quotients, so that 




-A^ T 
_j_ I j-i 

Sk -Sk 


A^ 

_ l_ 

Sk 


-A^ 


1 + 


+ \ 1 




= 0{ek), 


since 


1+ TT +W1+I 


< 1 and ——= 0{ek) by (33)-(34). Further, by Definition 


4.1 


O 


Cfc 


= o(i£,r). 


Thus, we can write 




CT£fe|Cfc| 


+ 0(|£,.r)=e. 


V 


£k 


+ W]sfe^ {qj) 


\ 




1 

i-f 

J 

£k 



We identify the right hand side as a product of a function w evaluated on the piecewise constant 
interpolation 7 ))'^ times a function z evaluated on the piecewise linear interpolation 7 ))^ that is, 

for all X G (x^ 1 , A* 1 ) there holds 
V 1-5 1 + 5 / 

(37) 


CTEfclCfcl 


= ei 


D. 


, (e'^ (a^,7r(a^))-P’(ay7f(a^))) +VAre^ (a:, 7 P"(a;)) (^ 7 P‘(a;) 


where is the centred difference quotient 


DeVix) := 


y{x + ^) - y{x - §) 


and 


m--= Pf{0 ■■= + 


i + ier 


We remark that the right-hand of (37) side is globally defined for x G {—£,£), and piecewise 
constant. Given an arbitrary fixed function y G Jig (—£, £), we can multiply (37) with y and 
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integrate by parts to obtain 

^Lj (sfe, 


(38) 

M-1 

E 

i=-M+i 


y{x) 


J-2 


I Cfc I 


■o(k,r) 


da; 


£i 

2 


y{x) 


D, 


da; 


(e^ (a;, iTix)) F (^ 7 ^( 2 :)) ) + Vatc^ (a;, 7r(a:)) Pf (^7^(a:) 
{Deuy{x)) {x, 7r(a;)) F (a 77 fc‘(a;)) + y{x)V nb^ (x, Pp ( 377 ^( 2 ^)) 


dx. 


We observe that the right-hand side does not depend on cr. Now recall that stoppage of the 
algorithm on level k implies that ALj{ek,CkO'ei) > 0 for both choices of the sign of cr = ± 1 . 

Let us assume for the moment that piecewise constant interpolation 7 ))'^ and the piecewise linear 
interpolation 7 ^' converge strongly to the same limit 7 as fc —>■ 00 in the sense that 


(39) 
and 

(40) 


hr-711^2^0, 


S7f - S7 


L 2 


0 . 


Then the argument can be finished as follows. Observe that ALj^Sk-, Ck^xBi) is non-negative for 
both choices of cr = ±1, whereas the last line of (38) does not depend on the chosen cr anymore. 




On each interval (^y_i j, we first choose cr to have the same sign as ^ y{x) dx. Then 

the sum on the left-hand side is non-negative. We pass to the limit fc —)■ c» in ([M|) and find 


M-l 


(41) 0 < lim 'V 

k—>-oo 


i=-M+l 


lx , ^ <xek\Ck\ 


= Bi 


[- (s^h)) h, 7 (a;)) F {^l{x)) + y{x)V nb^ {x,-i{x)) Pp (gy 7 (a;))] dx; 


here the convergence of the difference quotient Dg^,y{x) follows from [TT] Lemma 7.24]. 

Similarly, we then choose cr to have always the opposite sign and obtain the reversed inequality. 
Together this yields 


0 = Bi 


[- {-tyi^)) ( 2 ;, 7 ( 2 ;)) F (^ 7 ( 2 ;)) + y{x)VNB^ (x, 7 (x)) Pp (^ 7 ( 2 ;))] dx. 


As we can test all normal directions e^, with f = 2 ,..., n, in this fashion, we obtain the vectorial 
identity 


(42) 


0 = 


[- i^yi^)) (2^: 7 ( 2 ;)) F (^ 7 ( 2 :)) + y(x)Vjve'‘ (x, 7 (x)) Pp (^ 7 ( 2 :))] dx 


pn—1 


. By substituting the definitions of F and Pp, we recover the claim. 


□ 


4.3.1. Proof of (39) and (40). The previously assumed convergence of the interpolants is estab¬ 
lished in the following arguments. 

Lemma 4.6 (Estimates for interpolants). Given a function 7 '^ G FI^, let he the piecewise 
constant interpolation on an equidistant grid of size s, and similarly 7 P* the piecewise linear in¬ 
terpolation. Then the following error estimates hold 


- 7'^IL2 < 




L 2 
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and 


I A^pi, 
I da? ' 


da? ^ I 


l2 < Ce 


da?2 ' 


L 2 


Proof. This is a standard argument in Finite Elements, see for example [H Theorem 0.8.7 and 
Section 4]. □ 

For a fixed refinement level k, let qj be the point set associated to 7 ^, the output of the 
Birkhoff map at level k. Now, we first construct a quadratic interpolation of {q^} with 
the special property that its piecewise constant and piecewise linear interpolation coincide with 
the interpolations 7 ^'^ and 7 ^* introduced before. Specifically, in our situation the quadratic 


interpolation can be chosen in such a way that 


_di pq 
dx^lk 


L 2 


is bounded independently of k. 


Explicitly, to construct any two neighbouring nodes qj = and 

{X^ +£k,Y^ -\- dJf ) are connected via two quadratic splines s± on 1 ^ and 

with matching conditions 


X 




j+i’ J'+i 


S-{X^)=Y^ 


s'_{X'f) = 6 _ := ^ 


A’;, 


Set 


s+(X^+i) = X' 


k 

i+1’ 


sV(XjVi) = b+ := 


Ak 

^j+1 


■ A? 




s-{X^+i) = S+iX^rJ, s'_(Xf^j) = 4(Xf^i). 


j+ 


0+i 




On the boundary, we vary the definition slightly for j = —Mk and j = Mk — 1 respectively. 


.{x'im,) = y’: 


and 


s+(x^j = y, 


Mk' 


k 

Mfc? 


. A'? 

'_{XtMk) = b- :=-^ 


£fe 


>(XL) = b+ ■.= 


A 


Mk — 1 

Sk 


We recall that the polygons = ( 2 ^, fk{x)) obtained by the Birkhoff refinement satisfy the finite 
difference estimates (33)-(34). For the chosen assignments b± of the first derivatives at the nodes, 
the piecewise quadratic interpolation 7 ^'^ is continuously differentiable throughout. Further, it 
can be computed that the bound on the second difference quotients implies a bound 


(43) 


da:2 7fc 


<C, 


which holds uniformly in k. Ffence the sequence { 7 ^'^} is uniformly bounded in C^’^, and we infer 
convergence in to 7 G by Arzela-Ascoli. This implies that 

(44) hr-TlL^^O, 

and 


I A^pq - 

I da; fk 


da? ' I 


L 2 


0 


(45) 

as fc —>■ 00 . Choosing 7 ))'^ as the function 7 ^ in Lemma 4.6 we infer for all k that 


ii7r-7L.<ii7r-7riiE. + ii7r-7L. 


<Cel 


AA^pq 

da?2 !k 


L 2 


+ ii7r-^ 




and 


d pi 
dS7fc 


— A 
da? ^ 


L 2 


< iis7r 


-4-'vPq|| 

da? |Il2 


1=7’ 


—aII 

dx '11^2 


A Csk 


jB pq 

da?2 'k 


L 2 


I A-^pq. 

I dx Ik 


—aII 

da? MlL2 ■ 
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All expressions tend to 0 by (43)-(45). This proves the previously assumed claims (391 and (40). 


5. Numerical investigations 

Here we present experimental convergence rates and estimates for the computational effort for 
a simple benchmark problem. First, we obtain an explicit solution for a special metric which 
will serve as comparison for the numerical approximations on different discretisation levels £k- 


5.1. Special analytic solution. From the weak formulation of the geodesic equation for graphs (35), 
we obtain by integration by parts the following strong version 


ViveVl + l/f - + f ■ V^,eyi + |/f) ^ 


/' 


1/1 


^ _-13 J ; 

i + l/f 


which has to satisfied by a geodesic connecting the two points (±|,0) being a graph 'y{x) = 
(xJix)), /(±|) = 0. 

If we assume further that h does not depend on x, the last equation simplifies to 


0 = 




(Vjv/i (l + Iff) - /") . 


That is, we need to solve the simpler equation 

/" = (l + l/f) . 

Let us restrict ourselves to linear h, not depending on the ‘horizontal’ co-ordinate x, that is, 
h = h{y) = —arfy with a > 0 and n G a unit vector, |n| = 1. Then Vjyh = —an, hence, 

/" = —a -I- 1/1^^ n. With the ansatz f = fn, the equation can be rewritten as 

f'n = —an ^1 -I- \ ff) , 
hence we are left with the scalar equation 

f = -a (l + Iff) 

for f alone. Its solution is f{x) = tan(C' — ax). Thus, integrating once again we deduce 
fx) = ^lncos(—C -I- ax) + D. Matching the boundary conditions, we hnd that the geodesic 
connecting the two points (±|,0) is 

lix) = (x, fix)), a; e (-|, |), fix) = nfx) 

with 

, , 1 cos ax 

4>ix) = - In- j. 

a cosa| 

We remark that this representation requires a| < f, or equivalently a < j. Note that for 
a —>■ 0, the metric approaches the constant Euclidean metric h = 0, and we obtain the straight 
line segment (a;, / = 0) in the limit. 

For the other extreme, we remark that if a —)■ j, the geodesic converges to the set of two 
parallel lines {(±^,j/) | i/ > 0} connected at infinity. In fact, this set of constant and finite length 
will be the minimising configuration for all values of a > j. 
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Figure 3. Exponential metric, n = 2: fc = 1,..., 8. Left panel: The rate of the 
error proportional to e™. Right panel: The effort grows approximately with the 
rate 


5.2. Computational effort. We use the analytic solution found above as a benchmark test for 
convergence of the method. Let 

9ir{x,y) ■■= with h{y) = -ay. 

We now present computations in with £ = 2, so connecting (±1,0). Let us consider three 
different metrics by choosing three values oi = 0.65, 02 = 0.9, and 03 = 1.1, all less than |. 
Larger values of a emphasise the difference of g compared to the flat Euclidean metric, hence one 
expects a larger curvature in the geodesic, which bends further away from the horizontal straight 
connection between the boundary points (±1,0). We also explore how the performance of the 
algorithm depends on this geometrical feature. 

We begin by presenting numerical statistics of the geodesic computation. The plots in Figure]^ 
show the computational error of the polygonal approximation compared to the explicit geodesic 
(measured in the norm) in the left panel and the computational effort, measured in the number 
of affirmative Birkhoff steps, in the right panel. 

The error decreases as Sk = 2“^ decreases for fc = 1,... , 8 . In fact, the calculation shows the 
error is proportional to e™, where the exponent m varies from 0.40 to 0.53, depending on the 
choice of a in the metric. 

The effort increases as e decreases, and the simulation shows that the effort turns out to be 
(inversely) proportional to 1/e”^, where m ranges from 1.71 to 2.2, as a is taken from ai = 0.65, 
02 = 0.9, 03 = 1.1. 
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Figure 4. Exponential metric, n = 2: approximate solutions compared to the 
analytic solution (solid line). 


Hence, for geodesics of larger curvature variation, the Birkhoff procedure starts closer but 
converges with smaller rate in e, whereas the effort increases with largely uniform rate in e of 
about 1 /e^. 

Figure]^ shows how the polygons on the discrete level Sk approach the exact solution 


7 (x) 



COS ax 
cos a 


X e (-1,1). 


We show the approximation for the three different values for a simultaneously in one graph. In 
each case, the solid line is the exact geodesic, the dotted line is the polygonal approximation for 
63 = 1/8, with the thicker dots indicating the stencil points. The intermediate polygon with 
densely distributed stencils is the approximation for e^ = 1/128. 
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